ElectionForecast
load data
# description of data is here:
# http://www.stat.columbia.edu/~gelman/book/data/presidential.asc
electiondata <- read.table("election.txt", header = TRUE)visualizations
namescol <- names(electiondata)
# Fig 1: visualize response and covariates
for(j in 1:dim(electiondata)[2]){
hist(electiondata[, j], main = namescol[j], xlab = toString(round(summary(electiondata[, j]), 3)), ylab = "")
}
Remove NA’s
idxNA <- which(unlist(apply(electiondata, 1, function(x){sum(is.na(x))})) > 0)
print(electiondata[idxNA, ])## Dvote year state evotes constant n1 n2 n3 n4 s1 s2
## 251 0.260 1972 1 9 1 0.3711 -0.62 -0.62 -2.17 0.077 NA
## 301 0.573 1968 1 10 1 0.4189 0.49 0.00 1.30 NA 0.073
## 324 0.630 1968 24 7 1 0.4189 0.49 0.00 1.30 -0.484 NA
## 351 NA 1964 1 10 1 0.6915 0.74 0.74 1.12 0.073 0.167
## 352 0.659 1964 2 3 1 0.6915 0.74 0.74 1.12 -0.010 NA
## 361 0.788 1964 11 4 1 0.6915 0.74 0.74 1.12 -0.001 NA
## 374 0.129 1964 24 7 1 0.6915 0.74 0.74 1.12 NA 0.282
## 402 0.491 1960 2 3 1 0.4946 -0.68 0.00 -0.39 NA NA
## 411 0.500 1960 11 3 1 0.4946 -0.68 0.00 -0.39 NA NA
## 424 NA 1960 24 8 1 0.4946 -0.68 0.00 -0.39 0.282 0.158
## 451 0.589 1956 1 11 1 0.4409 -0.69 -0.69 -0.46 0.178 NA
## 453 0.534 1956 4 8 1 0.4409 -0.69 -0.69 -0.46 0.114 NA
## 458 0.427 1956 9 10 1 0.4409 -0.69 -0.69 -0.46 0.004 NA
## 459 0.666 1956 10 12 1 0.4409 -0.69 -0.69 -0.46 0.251 NA
## 466 0.426 1956 18 10 1 0.4409 -0.69 -0.69 -0.46 0.083 NA
## 472 0.704 1956 24 8 1 0.4409 -0.69 -0.69 -0.46 0.158 NA
## 481 0.507 1956 33 14 1 0.4409 -0.69 -0.69 -0.46 0.093 NA
## 488 0.643 1956 40 8 1 0.4409 -0.69 -0.69 -0.46 0.061 NA
## 490 0.497 1956 42 11 1 0.4409 -0.69 -0.69 -0.46 0.052 NA
## 491 0.443 1956 43 24 1 0.4409 -0.69 -0.69 -0.46 0.088 NA
## 494 0.409 1956 46 12 1 0.4409 -0.69 -0.69 -0.46 -0.011 NA
## 499 0.649 1952 1 11 1 0.4211 0.28 0.00 -0.24 NA 0.279
## 501 0.560 1952 4 8 1 0.4211 0.28 0.00 -0.24 NA 0.163
## 506 0.450 1952 9 10 1 0.4211 0.28 0.00 -0.24 NA 0.165
## 507 0.697 1952 10 12 1 0.4211 0.28 0.00 -0.24 NA 0.288
## 514 0.529 1952 18 10 1 0.4211 0.28 0.00 -0.24 NA 0.268
## 520 0.604 1952 24 8 1 0.4211 0.28 0.00 -0.24 NA 0.398
## 529 0.539 1952 33 14 1 0.4211 0.28 0.00 -0.24 NA 0.129
## 536 0.507 1952 40 8 1 0.4211 0.28 0.00 -0.24 NA 0.413
## 538 0.498 1952 42 11 1 0.4211 0.28 0.00 -0.24 NA 0.069
## 539 0.468 1952 43 24 1 0.4211 0.28 0.00 -0.24 NA 0.272
## 542 0.435 1952 46 12 1 0.4211 0.28 0.00 -0.24 NA 0.087
## 547 NA 1948 1 11 1 0.4561 0.39 0.39 1.78 0.279 0.306
## s3 s4 s5 s6 s7 s8 s9 r1 r2 r3
## 251 0 0 0.4811321 -1.23101099 -0.5169 -0.21445981 0.00000000 0 0 0
## 301 0 0 0.5000000 -1.08320756 -0.5169 -0.23250000 0.00000000 0 0 0
## 324 0 0 0.4836066 -0.08320756 -0.7610 -0.23250000 0.00000000 0 0 0
## 351 0 0 0.4811321 2.44829330 -0.5169 0.31750000 0.00000000 0 -1 -1
## 352 0 0 0.0000000 3.44829330 0.2181 0.14523025 0.00000000 0 0 0
## 361 0 0 0.2843137 0.84829330 0.5367 0.31750000 0.00000000 0 0 0
## 374 0 0 0.5000000 0.94829330 -0.7610 0.31750000 0.00000000 0 -1 -1
## 402 0 0 0.3500000 -2.00906169 0.2181 0.03773025 -0.06447399 0 0 0
## 411 0 0 0.1470588 -5.20906169 0.5367 0.14500000 0.10097422 0 0 0
## 424 0 0 0.5000000 3.29093831 -0.7610 0.14500000 -0.19114162 0 0 0
## 451 0 0 0.5000000 -4.75645805 -0.5169 -0.08458906 0.00000000 0 0 0
## 453 0 0 0.4797980 -2.65645805 -0.2121 0.13500000 0.00000000 0 0 0
## 458 0 0 0.4368421 -8.15645805 -0.3438 0.06126786 0.00000000 0 0 0
## 459 0 0 0.4853659 -1.35645805 -0.4633 0.13500000 0.00000000 0 0 0
## 466 0 0 0.5000000 -2.35645805 -0.3767 0.13500000 0.00000000 0 0 0
## 472 0 0 0.5000000 -2.55645805 -0.7610 0.13500000 0.00000000 0 0 0
## 481 0 0 0.4243697 0.04354195 -0.5367 0.13500000 0.00000000 0 0 0
## 488 0 0 0.5000000 3.44354195 -0.5476 0.13500000 0.00000000 0 0 0
## 490 0 1 0.3061224 0.24354195 -0.2737 0.13500000 0.00000000 0 0 0
## 491 -1 0 0.5000000 1.24354195 -0.2663 0.13500000 0.00000000 0 0 0
## 494 0 0 0.4400000 1.54354195 -0.6992 -0.02295759 0.00000000 0 0 0
## 499 0 1 0.4905660 0.21605450 -0.5169 -0.09708906 0.00000000 0 0 0
## 501 0 0 0.4800000 -6.08394550 -0.2121 0.15250000 0.00000000 0 0 0
## 506 0 0 0.4684211 1.51605450 -0.3438 0.04876786 0.00000000 0 0 0
## 507 0 0 0.4951220 3.81605450 -0.4633 0.15250000 0.00000000 0 0 0
## 514 0 0 0.5000000 1.81605450 -0.3767 0.15250000 0.00000000 0 0 0
## 520 0 0 0.5000000 0.31605450 -0.7610 0.15250000 0.00000000 0 0 0
## 529 0 0 0.4250000 -1.58394550 -0.5367 0.15250000 0.00000000 0 0 0
## 536 0 0 0.5000000 7.01605450 -0.5476 0.15250000 0.00000000 0 0 0
## 538 0 0 0.3080808 -3.88394550 -0.2737 0.13418020 0.00000000 0 0 0
## 539 -1 0 0.4933333 4.81605450 -0.2663 0.15250000 0.00000000 0 0 0
## 542 0 0 0.4285714 2.31605450 -0.6992 -0.03545759 0.00000000 0 0 0
## 547 0 0 0.4905660 -0.47588378 -0.5169 -0.14708906 0.00000000 0 0 0
## r4 r5 r6
## 251 0 0 0
## 301 0 0 0
## 324 0 0 0
## 351 0 0 0
## 352 0 0 0
## 361 0 0 0
## 374 0 0 0
## 402 0 0 0
## 411 0 0 0
## 424 0 0 0
## 451 0 0 0
## 453 0 0 0
## 458 0 0 0
## 459 0 0 0
## 466 0 0 0
## 472 0 0 0
## 481 0 0 0
## 488 0 0 0
## 490 0 0 0
## 491 0 0 0
## 494 0 0 0
## 499 0 0 0
## 501 0 0 0
## 506 0 0 0
## 507 0 0 0
## 514 0 0 0
## 520 0 0 0
## 529 0 0 0
## 536 0 0 0
## 538 0 0 0
## 539 0 0 0
## 542 0 0 0
## 547 0 0 0electiondata <- electiondata[setdiff(1:dim(electiondata)[1], idxNA), ]
# vote for the Democrats
Y <- electiondata$Dvote
# predictors
X <- electiondata[, 2:dim(electiondata)[2]]Pair wise correlations
library(ggplot2)
for(j in 2:dim(electiondata)[2]){
plot(electiondata[, j], electiondata$Dvote, xlab = names(electiondata)[j], ylab = "Dvote")
}
are previous years predictive of current year?
years_uniq <- sort(unique(electiondata$year))
for(j in 1:(length(years_uniq) - 1)){
prev <- electiondata$Dvote[electiondata$year == years_uniq[j]]
new <- electiondata$Dvote[electiondata$year == years_uniq[j+1]]
if(length(prev) == length(new)){
plot(prev, new, xlab = years_uniq[j], ylab = years_uniq[j+1], main = '')
}
}
First Model: simple linear regression
linearregressionresult <- lm(Y~as.matrix(X))
summary(linearregressionresult)##
## Call:
## lm(formula = Y ~ as.matrix(X))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.094012 -0.024325 -0.000339 0.023995 0.113614
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.001e+00 3.006e-01 -3.330 0.000927 ***
## as.matrix(X)year 6.588e-04 1.547e-04 4.258 2.44e-05 ***
## as.matrix(X)state 2.805e-05 1.092e-04 0.257 0.797284
## as.matrix(X)evotes 4.304e-04 1.791e-04 2.404 0.016566 *
## as.matrix(X)constant NA NA NA NA
## as.matrix(X)n1 3.784e-01 2.601e-02 14.548 < 2e-16 ***
## as.matrix(X)n2 -1.397e-02 5.782e-03 -2.416 0.016011 *
## as.matrix(X)n3 4.879e-02 7.866e-03 6.203 1.11e-09 ***
## as.matrix(X)n4 2.839e-02 1.857e-03 15.285 < 2e-16 ***
## as.matrix(X)s1 3.047e-01 3.582e-02 8.508 < 2e-16 ***
## as.matrix(X)s2 2.633e-01 2.878e-02 9.148 < 2e-16 ***
## as.matrix(X)s3 4.587e-02 8.119e-03 5.650 2.61e-08 ***
## as.matrix(X)s4 1.619e-02 7.873e-03 2.056 0.040290 *
## as.matrix(X)s5 3.826e-02 9.907e-03 3.862 0.000126 ***
## as.matrix(X)s6 1.181e-03 4.242e-04 2.783 0.005574 **
## as.matrix(X)s7 3.042e-02 5.588e-03 5.443 7.96e-08 ***
## as.matrix(X)s8 3.780e-02 1.758e-02 2.150 0.032026 *
## as.matrix(X)s9 1.226e-01 4.296e-02 2.854 0.004483 **
## as.matrix(X)r1 5.518e-02 7.695e-03 7.171 2.48e-12 ***
## as.matrix(X)r2 8.827e-02 1.624e-02 5.435 8.34e-08 ***
## as.matrix(X)r3 1.823e-01 2.585e-02 7.055 5.35e-12 ***
## as.matrix(X)r4 6.847e-02 1.217e-02 5.627 2.95e-08 ***
## as.matrix(X)r5 6.834e-02 1.277e-02 5.351 1.29e-07 ***
## as.matrix(X)r6 6.145e-02 9.229e-03 6.658 6.86e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03617 on 538 degrees of freedom
## Multiple R-squared: 0.8622, Adjusted R-squared: 0.8566
## F-statistic: 153 on 22 and 538 DF, p-value: < 2.2e-16plot(linearregressionresult$fitted.values, Y)
lines(quantile(Y, c(0.01, 0.99)), quantile(Y, c(0.01, 0.99)), col = "red")
# Rstan
require(rstanarm)## Loading required package: rstanarm## Loading required package: Rcpp## rstanarm (Version 2.18.2, packaged: 2018-11-08 22:19:38 UTC)## - Do not expect the default priors to remain the same in future rstanarm versions.## Thus, R scripts should specify priors explicitly, even if they are just the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())## - Plotting theme set to bayesplot::theme_default().postlinearmodel <- stan_glm(Dvote ~ .-1, data = electiondata,
family = gaussian(link = "identity"))##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000697 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.97 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 87.9119 seconds (Warm-up)
## Chain 1: 61.1866 seconds (Sampling)
## Chain 1: 149.099 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 8.3e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.83 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 91.21 seconds (Warm-up)
## Chain 2: 62.8687 seconds (Sampling)
## Chain 2: 154.079 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.77 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 70.4719 seconds (Warm-up)
## Chain 3: 61.5945 seconds (Sampling)
## Chain 3: 132.066 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7.5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 84.4976 seconds (Warm-up)
## Chain 4: 61.0738 seconds (Sampling)
## Chain 4: 145.571 seconds (Total)
## Chain 4:summary(postlinearmodel)##
## Model Info:
##
## function: stan_glm
## family: gaussian [identity]
## formula: Dvote ~ . - 1
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 4000 (posterior sample size)
## observations: 561
## predictors: 23
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## year 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## state 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## evotes 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## constant -0.4 0.2 -0.7 -0.5 -0.4 -0.3 0.0
## n1 0.4 0.0 0.3 0.4 0.4 0.4 0.4
## n2 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## n3 0.0 0.0 0.0 0.0 0.0 0.1 0.1
## n4 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s1 0.3 0.0 0.2 0.3 0.3 0.3 0.4
## s2 0.3 0.0 0.2 0.2 0.3 0.3 0.3
## s3 0.0 0.0 0.0 0.0 0.0 0.1 0.1
## s4 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s5 0.0 0.0 0.0 0.0 0.0 0.0 0.1
## s6 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s7 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s8 0.0 0.0 0.0 0.0 0.0 0.1 0.1
## s9 0.1 0.0 0.0 0.1 0.1 0.2 0.2
## r1 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r2 0.1 0.0 0.1 0.1 0.1 0.1 0.1
## r3 0.2 0.0 0.1 0.2 0.2 0.2 0.2
## r4 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r5 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r6 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## sigma 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## mean_PPD 0.5 0.0 0.5 0.5 0.5 0.5 0.5
## log-posterior 1039.2 3.6 1031.3 1036.9 1039.5 1041.7 1045.1
##
## Diagnostics:
## mcse Rhat n_eff
## year 0.0 1.0 2703
## state 0.0 1.0 7401
## evotes 0.0 1.0 5229
## constant 0.0 1.0 2794
## n1 0.0 1.0 3038
## n2 0.0 1.0 3375
## n3 0.0 1.0 3168
## n4 0.0 1.0 3564
## s1 0.0 1.0 3681
## s2 0.0 1.0 4944
## s3 0.0 1.0 5584
## s4 0.0 1.0 7240
## s5 0.0 1.0 3857
## s6 0.0 1.0 6616
## s7 0.0 1.0 4315
## s8 0.0 1.0 3575
## s9 0.0 1.0 5303
## r1 0.0 1.0 4501
## r2 0.0 1.0 3507
## r3 0.0 1.0 4089
## r4 0.0 1.0 5094
## r5 0.0 1.0 6088
## r6 0.0 1.0 5032
## sigma 0.0 1.0 5517
## mean_PPD 0.0 1.0 3688
## log-posterior 0.1 1.0 1388
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).posteriorsampleslm <- as.matrix(postlinearmodel)
coeflm <- as.numeric(linearregressionresult$coefficients)
coeflm[5] <- coeflm[1]
cbind(coef(postlinearmodel), coeflm[-1])## [,1] [,2]
## year 3.414303e-04 6.587634e-04
## state 1.851498e-05 2.805306e-05
## evotes 4.147651e-04 4.304138e-04
## constant -3.844254e-01 -1.001246e+00
## n1 3.941203e-01 3.784376e-01
## n2 -1.284045e-02 -1.397224e-02
## n3 4.525827e-02 4.878971e-02
## n4 2.686466e-02 2.838677e-02
## s1 2.959038e-01 3.047399e-01
## s2 2.538682e-01 2.632691e-01
## s3 4.491983e-02 4.587272e-02
## s4 1.644348e-02 1.618528e-02
## s5 4.077731e-02 3.825868e-02
## s6 1.140069e-03 1.180717e-03
## s7 2.999306e-02 3.041943e-02
## s8 4.626688e-02 3.779558e-02
## s9 1.298037e-01 1.226186e-01
## r1 5.538342e-02 5.518309e-02
## r2 9.152500e-02 8.826796e-02
## r3 1.806605e-01 1.823343e-01
## r4 6.566946e-02 6.847465e-02
## r5 6.737347e-02 6.833733e-02
## r6 6.156956e-02 6.144558e-02Question: is the regression model above valid?
linearregressionresult2 <- lm(Y~as.matrix(X)[, -(1:3)])
summary(linearregressionresult2)##
## Call:
## lm(formula = Y ~ as.matrix(X)[, -(1:3)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.101423 -0.024965 -0.000357 0.024100 0.115080
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2832418 0.0121572 23.298 < 2e-16 ***
## as.matrix(X)[, -(1:3)]constant NA NA NA NA
## as.matrix(X)[, -(1:3)]n1 0.4151813 0.0249152 16.664 < 2e-16 ***
## as.matrix(X)[, -(1:3)]n2 -0.0118705 0.0058524 -2.028 0.043020 *
## as.matrix(X)[, -(1:3)]n3 0.0415290 0.0077887 5.332 1.43e-07 ***
## as.matrix(X)[, -(1:3)]n4 0.0250209 0.0017138 14.600 < 2e-16 ***
## as.matrix(X)[, -(1:3)]s1 0.2956735 0.0362343 8.160 2.36e-15 ***
## as.matrix(X)[, -(1:3)]s2 0.2449690 0.0289011 8.476 < 2e-16 ***
## as.matrix(X)[, -(1:3)]s3 0.0400832 0.0080632 4.971 8.95e-07 ***
## as.matrix(X)[, -(1:3)]s4 0.0146022 0.0079497 1.837 0.066785 .
## as.matrix(X)[, -(1:3)]s5 0.0404499 0.0099340 4.072 5.36e-05 ***
## as.matrix(X)[, -(1:3)]s6 0.0011482 0.0004312 2.663 0.007981 **
## as.matrix(X)[, -(1:3)]s7 0.0298943 0.0056833 5.260 2.08e-07 ***
## as.matrix(X)[, -(1:3)]s8 0.0572498 0.0174137 3.288 0.001076 **
## as.matrix(X)[, -(1:3)]s9 0.1489682 0.0433676 3.435 0.000638 ***
## as.matrix(X)[, -(1:3)]r1 0.0573191 0.0077921 7.356 7.07e-13 ***
## as.matrix(X)[, -(1:3)]r2 0.0919283 0.0164507 5.588 3.64e-08 ***
## as.matrix(X)[, -(1:3)]r3 0.1859974 0.0262716 7.080 4.50e-12 ***
## as.matrix(X)[, -(1:3)]r4 0.0594448 0.0122382 4.857 1.56e-06 ***
## as.matrix(X)[, -(1:3)]r5 0.0678408 0.0129952 5.220 2.55e-07 ***
## as.matrix(X)[, -(1:3)]r6 0.0641587 0.0093380 6.871 1.76e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03684 on 541 degrees of freedom
## Multiple R-squared: 0.8562, Adjusted R-squared: 0.8512
## F-statistic: 169.6 on 19 and 541 DF, p-value: < 2.2e-16plot(linearregressionresult2$fitted.values, Y)
lines(quantile(Y, c(0.01, 0.99)), quantile(Y, c(0.01, 0.99)), col = "red")
require(rstanarm)
postlinearmodel2 <- stan_glm(Dvote ~ constant + n1+n2+n3+n4+s1+s2+s3+s4+s5+s6+s7+s8+s9+r1+r2+r3+r4+r5+r6, data = electiondata,
family = gaussian(link = "identity"))## Warning in center_x(x, sparse): Dropped empty interaction levels: constant##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000236 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.36 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.00273 seconds (Warm-up)
## Chain 1: 1.44319 seconds (Sampling)
## Chain 1: 3.44593 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 6.6e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.66 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.89684 seconds (Warm-up)
## Chain 2: 0.945545 seconds (Sampling)
## Chain 2: 2.84238 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.79661 seconds (Warm-up)
## Chain 3: 1.0029 seconds (Sampling)
## Chain 3: 2.7995 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'continuous' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 7.5e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.75 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.01573 seconds (Warm-up)
## Chain 4: 0.976628 seconds (Sampling)
## Chain 4: 2.99236 seconds (Total)
## Chain 4:summary(postlinearmodel2)##
## Model Info:
##
## function: stan_glm
## family: gaussian [identity]
## formula: Dvote ~ constant + n1 + n2 + n3 + n4 + s1 + s2 + s3 + s4 + s5 +
## s6 + s7 + s8 + s9 + r1 + r2 + r3 + r4 + r5 + r6
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 4000 (posterior sample size)
## observations: 561
## predictors: 21
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## (Intercept) 0.3 0.0 0.3 0.3 0.3 0.3 0.3
## n1 0.4 0.0 0.4 0.4 0.4 0.4 0.5
## n2 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## n3 0.0 0.0 0.0 0.0 0.0 0.0 0.1
## n4 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s1 0.3 0.0 0.2 0.3 0.3 0.3 0.4
## s2 0.2 0.0 0.2 0.2 0.2 0.3 0.3
## s3 0.0 0.0 0.0 0.0 0.0 0.0 0.1
## s4 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s5 0.0 0.0 0.0 0.0 0.0 0.0 0.1
## s6 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s7 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## s8 0.1 0.0 0.0 0.0 0.1 0.1 0.1
## s9 0.1 0.0 0.1 0.1 0.1 0.2 0.2
## r1 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r2 0.1 0.0 0.1 0.1 0.1 0.1 0.1
## r3 0.2 0.0 0.1 0.2 0.2 0.2 0.2
## r4 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r5 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## r6 0.1 0.0 0.0 0.1 0.1 0.1 0.1
## sigma 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## mean_PPD 0.5 0.0 0.5 0.5 0.5 0.5 0.5
## log-posterior 1035.2 3.2 1028.1 1033.3 1035.6 1037.6 1040.7
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 0.0 1.0 3360
## n1 0.0 1.0 3289
## n2 0.0 1.0 2839
## n3 0.0 1.0 3077
## n4 0.0 1.0 3085
## s1 0.0 1.0 3213
## s2 0.0 1.0 4103
## s3 0.0 1.0 6057
## s4 0.0 1.0 5568
## s5 0.0 1.0 3341
## s6 0.0 1.0 6510
## s7 0.0 1.0 3665
## s8 0.0 1.0 3117
## s9 0.0 1.0 4230
## r1 0.0 1.0 3820
## r2 0.0 1.0 3672
## r3 0.0 1.0 3533
## r4 0.0 1.0 3912
## r5 0.0 1.0 4726
## r6 0.0 1.0 3813
## sigma 0.0 1.0 5148
## mean_PPD 0.0 1.0 4992
## log-posterior 0.1 1.0 1672
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).posteriorsampleslm2 <- as.matrix(postlinearmodel2)
coeflm2 <- as.numeric(linearregressionresult2$coefficients)
coeflm2[2] <- coeflm2[1]
cbind(coef(postlinearmodel2)[1:20], coeflm2[-1])## [,1] [,2]
## (Intercept) 0.283152176 0.283241772
## n1 0.415314973 0.415181305
## n2 -0.012039292 -0.011870470
## n3 0.041399567 0.041529023
## n4 0.025079846 0.025020909
## s1 0.295345496 0.295673462
## s2 0.245190699 0.244968991
## s3 0.040091186 0.040083235
## s4 0.015208162 0.014602152
## s5 0.040263185 0.040449916
## s6 0.001139342 0.001148189
## s7 0.029693174 0.029894283
## s8 0.057265716 0.057249766
## s9 0.149750055 0.148968199
## r1 0.057097953 0.057319109
## r2 0.092486481 0.091928284
## r3 0.184502348 0.185997421
## r4 0.059437167 0.059444787
## r5 0.067463586 0.067840815
## r6 0.064119824 0.064158736How these two regression models differ?
# estimated sigma: residue s.d.
summary(posteriorsampleslm2[, 21])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03306 0.03613 0.03688 0.03690 0.03766 0.04136summary(posteriorsampleslm[, 24])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03254 0.03557 0.03630 0.03633 0.03705 0.04083Model checking
# posterior correlations of parameters
library(ggplot2)
library(reshape2)
ggplot(data = melt(cov2cor(vcov(postlinearmodel))), aes(x=Var1, y=Var2, fill=value)) +
geom_tile() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
# posterior predictive quantities
years <- sort(unique(electiondata$year))
ypred <- posteriorsampleslm[, 1:23] %*% t(X) + rnorm(4000) * posteriorsampleslm[, 24]
# examine residues
residuepred <- ypred - posteriorsampleslm[, 1:23] %*% t(X)
residueobs <- t(apply(posteriorsampleslm[, 1:23] %*% t(X), 1, function(z) Y-z))
par(mfrow = c(1, 2))
hist(residueobs, xlim = c(-0.15, 0.15), freq = FALSE)
hist(residuepred, xlim = c(-0.15, 0.15), freq = FALSE)
par(mfrow = c(1, 2))
for(j in 1:length(years)){
idx <- which(electiondata$year == years[j])
hist(apply(ypred[, idx], 1, mean), xlab="mean", freq=FALSE, main = years[j])
abline(v=mean(Y[idx]), col = "red")
hist(apply(ypred[, idx], 1, var), xlab="var", freq=FALSE, main = years[j])
abline(v=var(Y[idx]), col = "red")
}
# posterior correlations of parameters
ggplot(data = melt(cov2cor(vcov(postlinearmodel2))), aes(x=Var1, y=Var2, fill=value)) +
geom_tile() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
# posterior predictive quantities
years <- sort(unique(electiondata$year))
ypred2 <- posteriorsampleslm2[, 1:20] %*% t(X[4:23]) + rnorm(4000) * posteriorsampleslm2[, 21]
residuepred2 <- ypred2 - posteriorsampleslm2[, 1:20] %*% t(X[4:23])
residueobs2 <- t(apply(posteriorsampleslm2[, 1:20] %*% t(X[4:23]), 1, function(z) Y-z))
par(mfrow = c(1, 2))
hist(residueobs2, xlim = c(-0.15, 0.15), freq = FALSE)
hist(residuepred2, xlim = c(-0.15, 0.15), freq = FALSE)
par(mfrow = c(1, 2))
for(j in 1:length(years)){
idx <- which(electiondata$year == years[j])
hist(apply(ypred2[, idx], 1, mean), xlab="mean", freq=FALSE, main = years[j])
abline(v=mean(Y[idx]), col = "red")
hist(apply(ypred2[, idx], 1, var), xlab="variance", freq=FALSE, main = years[j])
abline(v=var(Y[idx]), col = "red")
}
Year-to-Year, State-to-State
# does last year predict this year?
par(mfrow = c(2, 2))
for(k in 2:length(years)){
idxcurrent <- which(electiondata$year == years[k])
idxpast <- which(electiondata$year == years[k-1])
commonstates <- intersect(electiondata$state[idxcurrent], electiondata$state[idxpast])
idxcurrent <- idxcurrent[which(electiondata$state[idxcurrent] %in% commonstates)]
idxpast <- idxpast[which(electiondata$state[idxpast] %in% commonstates)]
plot(Y[idxpast], Y[idxcurrent], ylab = years[k], xlab = years[k-1])
}
# does other states predict current state?
states <- sort(unique(electiondata$state))
Rmat <- array(NA, c(length(states), length(years), 2))
for(k in 1:length(states)){
for(j in 1:length(years)){
idxcs <- intersect(which(electiondata$year == years[j]), which(electiondata$state == states[k]))
idxos <- setdiff(which(electiondata$year == years[j]), idxcs)
if(length(idxcs) > 0 && length(idxos) > 0){
Rmat[k, j, ] <- c(mean(Y[idxcs]), mean(Y[idxos]))
}
}
}
par(mfrow = c(3,3))
for(i in 1:length(states)){
plot(Rmat[i, , 1], Rmat[i, , 2])
}
# Expand Model
# a varying coefficient model
Ymat <- array(0, c(length(states), length(years)))
Xmat <- array(0, c(length(states), length(years), 14))
for(i in 1:length(years)){
idx <- which(electiondata$year == years[i])
Ymat[electiondata$state[idx], i] <- Y[idx]
for(j in 1:14){
Xmat[electiondata$state[idx], i, j] <- electiondata[idx,j+4]
}
}
Southlabel <- rep(0, length(states))
SouthStates <- c(1, 4, 9, 10, 17, 18, 24, 33, 36, 40, 42, 43, 46)
Southlabel[SouthStates] <- 1
InputDataRegion <- list(Nt = length(years), Ns = length(states), p = 14, Y = Ymat, X = Xmat, Southlabel = Southlabel)
require(rstan)## Loading required package: rstan## Loading required package: StanHeaders## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)samplesRegionalModel <- stan(file = "RegionElection.stan", data = InputDataRegion)##
## SAMPLING FOR MODEL 'RegionElection' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000318 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.18 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 133.834 seconds (Warm-up)
## Chain 1: 158.751 seconds (Sampling)
## Chain 1: 292.585 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'RegionElection' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 0.000144 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.44 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 133.466 seconds (Warm-up)
## Chain 2: 158.844 seconds (Sampling)
## Chain 2: 292.31 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'RegionElection' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 0.000205 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.05 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 138.68 seconds (Warm-up)
## Chain 3: 159.796 seconds (Sampling)
## Chain 3: 298.476 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'RegionElection' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 0.000155 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.55 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 143.895 seconds (Warm-up)
## Chain 4: 165.493 seconds (Sampling)
## Chain 4: 309.388 seconds (Total)
## Chain 4:## Warning: There were 4000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problemssamplesRM <- extract(samplesRegionalModel)Posterior check
whichgammas <- Southlabel[electiondata$state] + 1
gammasamples <- sapply(1:length(whichgammas), function(j){
samplesRM$gamma[, whichgammas[j], which(electiondata$year[j] == years)]
})
deltasamples <- sapply(1:length(whichgammas), function(j){
samplesRM$delta[, which(electiondata$year[j] == years)]
})
Xbetasamples <- samplesRM$beta %*% t(electiondata[,5:18])
# calculate fitted values
fittedvaluesRM <- Xbetasamples + gammasamples + deltasamples
# how is the fit?
plot(fittedvaluesRM[1, ], Y)
# posterior predictive quantities
predictiedvaluesRM <- fittedvaluesRM + t(sapply(samplesRM$sigma, function(s) rnorm(length(whichgammas)) * s))
# residuals
residuepredRM <- predictiedvaluesRM - fittedvaluesRM
residueobsRM <- t(apply(fittedvaluesRM, 1, function(x) Y - x))
par(mfrow = c(1, 2))
hist(residueobsRM)
hist(residuepredRM)
# predictive means and variances VS observed
par(mfrow = c(1, 2))
for(j in 1:length(years)){
idx <- which(electiondata$year == years[j])
hist(apply(predictiedvaluesRM[, idx], 1, mean), xlab="mean", freq=FALSE, main = years[j])
abline(v=mean(Y[idx]), col = "red")
hist(apply(predictiedvaluesRM[, idx], 1, var), xlab="variance", freq=FALSE, main = years[j])
abline(v=var(Y[idx]), col = "red")
}
Check robustness to model specification
# Robustness to priors
# try other priors
# Robustness to Gaussian assumption
# try student t distributionModel selection and model averaging
# Are all the predictors useful?
# What is the probability that predictor j is not useful, i.e. beta_j = 0?